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Abstract 



00 

^vq We present a new, model-independent approach for measuring non-Gaussianity of the Cosmic 

Microwave Background (CMB) anisotropy pattern. Our approach is based on the empirical dis- 
tribution function of the normalized spherical harmonic expansion coefficients a£ m of a nearly 

O ' 

full-sky CMB map, like the ones expected from forthcoming satellite experiments. Using a set of 

\Q _ 

Kolmogorov-Smirnov type tests, we check for Gaussianity and independency of the at m . We test 
the method on two non-Gaussian toy-models of the CMB, one generated in spherical harmonic 
Qlt space and one in pixel (real) space. We also provide some rigorous results, possibly of independent 

2 ■ interest, on the exact distribution of the spherical harmonic coefficients normalized by an estimated 

C$ , angular power spectrum. 
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I. INTRODUCTION 



Temperature fluctuations in the CMB are an invaluable tool to constraint cosmological 
models and the process of structure formation in the universe. According to the standard 
version of the most popular theory of the very early universe, the so-called theory of 
inflation, the density fluctuations at primordial epochs should be Gaussian or very close to 
Gaussian distributed (see the reviews 0,0). This implies that the temperature fluctuations 
in the CMB as observed today should also be close to Gaussian, because they are related by 
linear theory to fluctuations in the early universe. However, the details of the inflationary 
scenario are still rather unclear. Newer, more sophisticated, versions of inflation predict 
small deviations from Gaussianity @, |], [5], |j, 0, ||. Detecting this non-Gaussianity in the 
CMB would thus be important for understanding the physics of the very early universe. 
There are also other possible mechanisms for the creation of non-Gaussianity in the CMB. 
If the universe has undergone a phase transition at early times, this could have given rise to 
topological defects (See Ref. |J for a review). These defects would show up as non-Gaussian 
features in the CMB temperature fluctuation field. 



With the high resolution data from the satellite missions MAP|33| and Planck 



Surveyor ||34|| one could be able to detect possible deviations from non-Gaussianity in the 
CMB. If this is indeed detected it would have a big impact on our understanding of the 
physics of the early universe. For these reasons, the search for procedures to test for 
Gaussianity in high resolution data has recently drawn an enormous amount of attention 
in the CMB literature. A number of methods were proposed, many based upon topological 
properties of spherical Gaussian fields: the behavior of Minkowski functionals [|H], 111] , 



temperature correlation functions |L2], the peak to peak correlation function p3|| , skewness 
and kurtosis of the temperature field |14| and local curvature properties of Gaussian and 
chi-squared fields Other works have focussed on harmonic space approaches: analysis 



of the bispectrum and its normalized version [l7j and the bispectrum in the flat sky 
approximation []I8| . The explicit form of the trispectrum for CMB data was derived in 
|19|, pO|j. Applications to COBE, Maxima and Boomerang data have also drawn enormous 



attention and raised wide debate EIL E2L E3l E4 
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Another reason to look for non-Gaussianities in CMB data, is to detect effects of 
systematic errors in the CMB map. One of these systematic effects can be stripes from the 
1/f noise in the detectors which have not been properly removed. Another systematics 
which could induce non-physical non-Gaussianities in the map could be the effect of 
straylight contamination from the galaxy or distortions of the main beam caused by the 
optics of the telescope. Finally astrophysical foregrounds like synchrotron emission, thermal 
dust emission or free-free emission from the galaxy and other extra galactic sources could be 
wrongly interpreted as a physical non-Gaussian signal. It is important that one uses data 
from simulated experiments to check whether these systematic effects induce non-physical 
non-Gaussian signatures in the CMB map. 



Our purpose here is to propose a new procedure to detect non-Gaussianity in harmonic 
space. More precisely, let T(9, if) denote the CMB fluctuations field, which we assume, 
as always, to be homogeneous and isotropic, for < # < 7r, < tp < 2n. Assuming that 
T(9, <p) has zero mean and finite second moments, it is well-known that the following spectral 
representation holds: 

oo I 

T(9, if) = 22 ^2 a tmYi m (9, y?) , 

£=l m=-l 

where Ye m (9,ip) denotes the spherical harmonics. The random coefficients (amplitudes) 
{aim} have zero- mean with variance (|a^ m | 2 ) = Cg. They are uncorrelated over t and \m\: 
(ag m a^, ml ) = Cg5f 8™ , and a^ m = a} _ m . The sequence {Cg} denotes the angular power 
spectrum of the random field and the asterisk complex conjugation. Furthermore, if T(9, (p) 
is Gaussian, the {ag m } have a complex Gaussian distribution. Upon observing T(9, <p) on 
the full sky, the random coefficients can be obtained through the inversion formula 



0>en 



/7T /*7T 
/ T (9, <p)Y; m (9,<p) sin 9d9d(p , m = 0,±1,...,±Z , / = 1,2,... . (1) 
-7T JO 



Our purpose is to study the empirical distribution function for the {a^ m } , and to use these 
results to implement tests for non-Gaussianity in harmonic space. We shall assume that 
the angular power spectrum is unknown, and the sequence {Ce} estimated from the data; 
as we show below, this has a nonnegligible effect on the behavior of the test, no matter 
how good is the resolution of the experiment. The plan of the paper is as follows: the 
procedure we advocate is described in Sections |D] and |T|; Section |V| presents some empirical 



results, whereas Section |V| is devoted to discussion and to directions for further work. Some 
mathematical results are collected in the Appendix. 



II. THE UNIVARIATE EMPIRICAL PROCESS 



To motivate our procedure, we start from the (unrealistic) assumption that the sequence 
of coefficients in the angular power spectrum of T(9,(p), i.e. {Ci} £=12 L , is known; here, 
L is the highest observable multipole, which depends upon the resolution and noise level 
of the experiment (for instance, L ~ 2000 for Planck). Now recall that, under Gaussianity, 
\ a m\ 2 /Ci and 2\a^ m \ 2 /C^ are standard chi-square variates, with one and two degrees of free- 
dom, respectively; furthermore, they are independent and identically distributed (i.i.d.) over 
all I. It is clearly computationally convenient to introduce a transformation, in order to work 
with random variables that have an uniform distribution on [0, 1]. This can be achieved as 
follows. Write $j(x), i — 1,2, for the cumulative distribution function of a chi-square with 
% degrees of freedom; recall that 

Jo V2nu 



$i(x) 



exp(—u/2)du 



f x 1 

,{x) = / - exp(— u/2)du — 1 - exp(— x/2) , <& 2 (<*) = — 2 l°g(l — a ) 
Jo 2 



Now introduce the Smirnov transformation (see for instance Ref . Pq] ) 

u eo = $i (&jp\ , u tm = $ 2 {^§Pj , m = 1, 2, 1,1 = 1, 2, ...L . (2) 

It is immediate to see that the random variables of the triangular array {ii£ m } are i.i.d. with 
a uniform distribution in [0, 1], i.e. 



P [u £m <a] = P 
= P 



$ 2 (^f- I < a 

Up 



$2 [^ 2 1 ( Qt )] = a 



and likewise for m^ , because $i, $2 are strictly increasing and hence invertible. For < a < 
1, we can hence define the empirical distribution function 



F f (a) 



1 - 



< a) 



m=0 
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l(x) denoting the indicator function, i.e the function that takes value unity if ue m smaller or 
equal than a, zero otherwise; Fe(a) evaluates the proportion of observed ug m which is below a 
certain value, and thus provides a sample analogue of P (u£ m < a) . Then, if the distribution 
of the U£ m is uniform indeed, that is, if the a^ m are actually Gaussian, the Glivenko-Cantelli 
Theorem [[25], |2Bj ensures that, with probability one, Fi(a) — > a, uniformly in [0,1]. This 



can be viewed as a most intuitive conclusion, as it states that the relative frequency with 
which a certain event occurs converges to the probability of the same event, as the number 
of experiments grows. It is therefore natural to look at the distance between Fg(a) and a for 
a test of the assumption that the ag m are Gaussian. More precisely, a centered and rescaled 
version of Fi(a) provides the (sequence of) empirical processes 

G t {ai) = V£ + l{F e (a) -a}, a G [0, 1] . 

The idea is that, if the ae m are not Gaussian, then Fe(a) —>■ F'(a), for some F'(a) such that 
sup \F'(a) — a\ > 0, and hence Gi{a) will exhibit very high values as a distinctive feature of 
non-Gaussianity, at least over some a. On the other hand, if the ae m are Gaussian, i^(ot) — a 
converges to zero for all a, whereas Gz{a) converges to a well-known process, the Brownian 



bridge |2^, [271 , the sup of which has a tabulated distribution which can be used to derive 
threshold values. 

An apparent deviation from Gaussianity over some £, however, must be careful assessed; 
indeed, in the Gaussian case, the single Gf (x) are independently distributed; hence, if we run 
a test at the 95% confidence level, we should expect approximately 5% of the multipoles 
to provide results above the threshold level. To avoid spurious detections, our aim is to 
combine rigorously the information over all different multipoles I into a single statistic. We 
shall hence focus on the partial sum process 

l [Lr] 

K L (r,a) = —Y^G e (a) , a G [0, 1] , r G [0, 1] , 

which has not been considered so far in the literature. The idea is to evaluate what random 
fluctuations we can expect over the multipoles, if the underlying field is Gaussian. Note 
that, if there is any departure from Gaussianity, then the behavior of Ki(r,a) should also 
detect the location of the non-Gaussianity in harmonic space. 
It is immediate to see that 

(K L (r,a)) = 0, 
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whereas simple calculations show also that, 

\im (K L (r lj a l )K L (r 2 ,a 2 )) — min(rx, r 2 ) {min(o;i, a 2 )(l - max(o;i, a 2 ))} ■ (3) 

L^oo 

With more effort it is actually possible to establish a much stronger result, i.e. functional 
convergence in the distribution; more precisely fl2~8|, as L — > oo 



K L (r,a) K(r,a) , 

where K(r, a) (the so-called Kiefer-Muller process) is a Gaussian zero mean random field 
on [0,1] x [0,1] (a and r in [0,1]) with covariance function given in (|3|). Here, =>- means 



weak convergence in a functional sense | 26| ; this is a much stronger notion of convergence 
than simply requiring that the distribution of K^r^a) is asymptotically the same as the 
distribution of K(r, a), for all fixed r and a. The difference can be explained as follows: for 
statistical inference, we must actually consider some functional of Ki(r, a), for instance the 
sup in a Kolmogorov-Smirnov type test. Now convergence in the distribution for every fixed 
point does not entail convergence of continuous functional such as the sup, whereas this is 
granted by the stronger notion of convergence we are exploiting here. 

We are now in a position to relax the restrictive assumption that the angular power spectrum 
is known, to consider the more realistic case where the latter is estimated from the data. 
We take 



Ce ~ 27TT ^ K 12 



lm | j 



so that a Smirnov-type transformation analogous to (0) now reads 



u eo = $i I l^-] , u lm = $ 2 ( 2 ^jH ) , m = 0, 1, 2, I, I = 1, 2. 



Cg / V C L 

Likewise, we have an estimated empirical distribution function 



1 £ 

and the empirical process with estimated parameters 

Gtia) = y/i+1 ^F t (a) -a} 
1 [Lr] ^ 

K L (a, r) = —= V G e (a) ,0<a<l,0<r<l. 
V L 



It is important to stress that, due to the presence of estimated parameters, the normalized 
random coefficients ug m are no longer independent as m varies, for a fixed l\ on the other 
hand, independence across different multipoles i is maintained. 

Because of course Ci converges to Cg as i grows, it might be conjectured that the effect 
of using estimated parameters becomes asymptotically negligible, for L large. In fact, we 
can show that this is not the case. It is shown in the Appendix that we have, as I — > oo 

( Uu im <a)) = a + ^- + o(]) , (4) 



( Uum<a))=a + o(jj , (5) 

where 

b(a) = (I - a) log(l - a) + hi - a) log 2 (l - a) . (6) 
Here o Q) means that the remaining terms go to zero faster than (|J. Also, as L — > oo, 

lim (K L (a,r)) = 2^/rb{a) . (7) 

L^oo 

Equation (^) shows how failing to take into account the presence of estimated parameters 
may results in unwarranted conclusions: the limiting field entails a non-vanishing bias and 
hence needs further centering before reliable inference can take place. In view of we 
have easily 

1 , ,L 



(G e (a)) = (F e (a)) - a \~ ^==b{a) + o( 

so that the bias is negligible for I large if one focuses on a single row, whereas considering the 
whole array, i.e. summing from £ = 1 to L, makes the bias nonnegligible in the aggregate. 

The presence of estimated parameters has also a nonnegligible effects on the behavior 
of the covariances, more precisely, it can be shown that we have, as L — > oo, and for all 
< a h a 2 < 1, < n,r 2 < 1, 

lim Cov \ K L (ai,ri), K L (a 2 , r 2 ) \ 
= min(ri, r 2 ) min(ai, a 2 ) {1 — max(ai, 0-2)} 

-min(ri,r 2 )[(f - «i)(I - a 2 ) log(l - ai) log(l - a 2 )] , (8) 



see Ref. 28 for further details. 
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In view of the previous results, for statistical inference we suggest to focus on the mean- 
corrected field 



K' L (a, r) = K L (a, r) - 2^Fb(a) . (9) 

An alternative strategy would be to implement the bias correction directly on Gg(a), by 
using the exact (rather than asymptotic) expression for the bias which is provided in the 
appendix. Our experience with simulated maps suggests that the difference between these 
two approaches is quite negligible as a first approximation. As before, it is possible to show 
a much stronger result than convergence of the mean and variance, i.e., as L — > oo 

K' L (a,r) => K(a,r) , 

where we write K(a,r) for the zero mean Gaussian process on [0, 1] x [0, 1] with zero mean 
and covariance given in (|8|); the convergence is meant in the uniform, functional sense 



introduced before |26|, |27fl. It is important to notice that the limiting distribution does not 
depend on any nuisance parameter, i.e. the asymptotic distribution is completely model- 
independent, given Gaussianity. The structure of the limiting covariance measure may have 
some interest for statistical application: for instance, it is well-known from the theory of 
Gaussian fields that maxima will approximately occur in the region where the variance 



(equations | and §) of the field itself peaks p9| . In our case, this corresponds roughly to 
a ~ 0.4, r ~ 1: a maximum located far away from this value may by itself suggest some 
non- Gaussianity in the a^ m . 

The previous results, however, afford much more powerful and well-established tests for 
Gaussianity. For instance, a Kolmogorov-Smirnov type of test is implemented, for any 
suitably large L, if we evaluate 

= sup sup \K' L (a,r)\ , (10) 

0<r<l 0<a<l 

and compare the observed value with the threshold level obtained, for any desired size of 
the test, by 

= sup sup \K'(a,r)\ ; 

0<r<l 0<a<l 

the latter value can be readily derived by Monte Carlo simulation; note that the limiting 
distribution does not entail any unknown, nuisance parameter. Likewise, a Cramer- Von 
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Mises type test can be implemented by looking at (see Ref. ||25||) 




Jo Jo 



The relative performance of ( |T0D and (|TT| ) depends on the nature of departure from Gaus- 
sianity; for instance, (|H]) may perform better in cases when there is a relatively strong 
non-Gaussianity, concentrated on a limited subsets of multipoles, whereas (|TT| ) is better 



suited for circumstances where the non-Gaussian behavior is more evenly spread over many 
different angular frequencies. Likewise many other types of goodness-of-fit statistics can be 
easily implemented, based upon the notion that the field K' L (a, r) should diverge, at least for 
some (a, r), if non-Gaussianity is truly present in the marginal distribution of the spherical 
harmonic coefficients. 

III. THE MULTIVARIATE EMPIRICAL PROCESS 

The procedure of the previous Section is centered upon the search for non-Gaussianity 
in the marginal distribution of the {(%«}. The empirical results presented in Section [IV] 
for some toy models show that the resulting tests are very efficient indeed, if the {ae m } are 
actually non-Gaussian. For some other toy models, namely those where the non-Gaussianity 
is strong in pixel space, the proposed procedure turns out to have much less power. This 
can be intuitively explained as follows. Take T(6, tp) to be, for instance, some nonlinear 
transform of an underlying Gaussian field, and evaluate the spherical harmonic coefficients 
by ([!]); then, due to a central limit theorem argument, the resulting linear combinations 
over the pixels can have a distribution much closer to Gaussian in harmonic space. In this 
case the non-Gaussianity may show up as a dependency between a^ m of different multipoles. 
With this in mind, in the present section we develop some more powerful tests for non- 
Gaussianity, which consider the joint distribution of the ae m coefficients over different i and 
m values. More precisely, in the bivariate case we look at the joint empirical distribution 
function 




mi=0 




9 



for some integer A m > 0; in the sequel, we shall always use the sum modulus £ 2 , e.g., 



m + A m = m + A r 



£2 

general multivariate case we focus on 



2, if m + A m > £ 2 (here [x] means integer- value) . In the 



1 h f k } 

F tl ...t h (ai, ...,a k ) = - ^ |l(^ im < ai) J|l(^ i)m+ A roi < , A mi > . (12) 

Again if the a^ m are Gaussian, F^^ax, oi 2 ) — > a\a 2 (or F £l ...e k (ai, ctk) — ► Yli=x a «) anc ^ 
we obtain the multivariate empirical process in a similar way as for the univariate process 

G il e 2 (a 1 ,a 2 ) = a/ {l x + 1) ^F hh (ax, a 2 ) - ai^j , 

and 

r * 

G^...4(ai, «fc) = V (4 + 1) < ^...^(ai, CKfc) - J ] 



Oii 
i=l 



In the sequel, it is convenient to write £x = £ and £i + \ = £ + A«, for some positive and 
distinct integers < A^, A^ 2 , A^-i << L. For < ai,a! 2 < 1, < r < 1, we hence 
focus on 

1 I( L - A «)' r l _ 
i^L(«i,«2,r) = —=== V G^ + A n (ai,a 2 ) , 

As discussed in Section || using the estimated gives a bias which needs to be corrected 
for. Some simple computations show that 

k k , 1 1 s 

(IlK^m, < ao> = n { a * + + °y J 

= n^ + 7S n a i 6 K)+oy. (13) 

i=l i=l j=l,jj=i 

This gives the bias corrected K' L (ax, 02, r) 

a 2 , r) = -K"l(«i, a 2 , r) - 2\/rai&(a! 2 ) - 2 v / ra 2 6(ai) . (14) 

For the general multivariate case, we take the integers A& to be strictly increasing, < 
A^i < Ai2 < ... < Af^k-i « L, and we consider 

1 [(£-A fc _0r] ^ 

A~l(q!i, ...,a k ,r) = — = G^.. /+Aek _ 1 (a l , a fc ) , 

V ~~ A^fc-i 



=1 
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k k 

K' l (olx, a k , r) = K L {a x , a k , r) - 2^) j j J ajb(ai) . (15) 

i=l j=l,j^i 

The rationale for these procedures can be explained as follows; although, as motivated before, 
the marginal distribution of the spherical harmonic coefficients can be close to Gaussian even 
in circumstances where T(8, if) is not, the joint assumption that the at m are Gaussian and 
independent uniquely identifies a Gaussian field in pixel space, i.e., it provides necessary 
and sufficient conditions. Therefore by looking at multiple rows a non-Gaussian feature can 
be more likely detected. 

By an argument similar to Section f|, if is possible to establish the weak convergence 
of K' L (aci, CKfc, r) to zero mean Gaussian fields with a complicated covariance function, 
which is not reported for brevity's sake. We stress that the limiting distribution changes 
when we vary the number of rows considered; on the other hand, for L large the effect 
of the spacing factors A^, A m is minimal (in the Gaussian case): hence, at least as a first 
approximation, a single Monte Carlo tabulation is sufficient, given the number of rows 
included. The robustness of the limiting result for varying A's is further investigated in the 
empirical section. 

In the previous discussion, we focussed for simplicity on the case where only a single 
spherical harmonic coefficient was selected from any multipole t. This assumption may be 
relaxed to consider the joint distribution of the spherical harmonic coefficients within the 
same row £. The limiting form of the bias term is here more complicated, however, due to 
the dependence among the normalized coefficients ui m over the same row £. To evaluate this 
bias, we can provide the following general result: for any < a%, a p < 1 such that 

>A 21og(l - Oj) 

^ 2£ + l ~ ' 

i=l 

(a condition which is always fulfilled, provided I is large enough), we have 

(flU« tmi < a,)) = (l + E 2 '°2f + T J> )' 1/2 ' (16) 

where T p is the class of all 2 P subsets 7 of {1, 2, p} , and #(7) denotes the number of 
elements in 7. For instance in the bivariate case p = 2 we have T p = {0, (1), (2), (1, 2)} and 
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we obtain, for t large enough, 



1 + 21og(l- ai )^ 1/2 ' oWi_^X^ 



2£+ 1 



/ 21og(l-a 2 ) Y 
V 2£+l J 



21og(l- Q!l ) 21og(l-« 2 )^ 
1 2£ + 1 2£ + 1 

111 I 

= ai"2 + -^"2& («i) + -«i6(a 2 ) + ^c(«i, a 2 ) + o(-) , 

where we defined 

c(ai, a 2 ) = - log (1 - Qfi) log (1 - a 2 ) (1 - «i)(l - a 2 ) • 
Likewise, for I large enough we have also 

( l(uimi - ai)l(uim 2 < a 2 )l(u £m3 < a 3 )) 

8=1 V 7 i=l \ j=l,j^i 

/ ^ 21o g (l- aj 

I ^ 2£+l 
\ j=i 

= aia 2 a3 + ^jOL X a 2 b (a: 3 ) + ^aia 3 6(a 2 ) + ja 2 a 3 b(a 3 ) 

1 1 1 

+-«ic(a 2 , a 3 ) + ^a 2 c(a 2 , a 3 ) + -a 3 c(a 2 , a 3 ) . 

In higher dimensions, the bias structure is readily provided by analogous computations, and 
we do not report it here for brevity's sake. 

From the previous results, it is simple to derive the exact distribution of the normalized 
(and squared) spherical harmonic coefficients, in the presence of estimated parameters. We 
note first that 

p 

(Y\_l(u£ mi < Qfi)) = P(u imi < Of!, ...,Ui mp < OCp) , 
i=l 

always, by the definition of the indicator function. Also, \a£ m \ 2 /Ce < 1 + ~, by construction. 
Hence, for < < l + |, we obtain from ([16|) 

p( K^P < KJ_ < } = y ( _ 1)#(7) / j _ y 



12 



For instance, in the bivariate case, for < x\, x-i < I + ~ we have 



2 

Lit L/c 



\ £-1/2 / o \ <-l/2 / o o \ ^-V2 



1- 1-^^ - 1 



2£+lJ V 2£+iy V 2 ^+ 1 2£+l j 

The previous expressions may have some independent interest when considering exact sta- 
tistical inference for CMB data (for instance, in likelihood analysis). 

IV. EMPIRICAL SECTION 

In this section we will test our method for two classes of toy models. First we will use 
a non-Gaussianity generated in spherical harmonic space which should be easily detected 
as our method is based on the spherical harmonic coefficients. In the second class of test 
models we will generate the non-Gaussianity in pixel space . For these models we expect a 
drop in the detection level at least for the univariate test, as the transformation to spherical 
harmonic space may make the spherical harmonic coefficients more Gaussian using a central 
limit theorem argument. 

To generate a sequence of non-Gaussian ci£ m , we consider first the same model as in 
,i.e. 



Ref. 30 



7 ng — qR _i_ 



where af m , a\ m are independent xl variables, normalized to have zero mean and the same 
angular power spectrum as the Gaussian part of the model. 
Our field is then defined as (Model 1) 

oo I 

e=i m =-l 

a lm = (l-/3)a £ G m + /3<f. 

Note that {a^} and {a^f} are independent and have by construction an identical angular 
power spectrum; the percentage of non-Gaussianity in the model is then uniquely determined 

as 

(3 2 

fNG = (1 -/?)' + ■ (17) 
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FIG. 1: The K max (univariate) from a Monte Carlo simulation of 500 Gaussian CMB realizations. 
The colored areas indicate the K mSuX intervals where 68%, 95% and 99% of the Gaussian realizations 
are found. If K mSuX is found in one of these intervals, this is taken as a 1, 2 and 3 sigma detection 
of non-Gaussianity respectively. 



Needless to say, this toy model is unphysical, but it is helpful to illustrate some characteristics 
of our approach. Note the bispectrum B™j^ m3 = (a£ imi ai 2m2 ae 3m3 ) here is identically zero, 
unless mi — m 2 — m 3 — 0; hence the selection rules implied by Wigner's 3j coefficients 
are not satisfied and the model is (slightly) anisotropic. We do not view this as a major 
difficulty, however, as the effect is minimal and this model is introduced here for purely 
expository purposes. 

We consider Komogorov-Smirnov type tests as discussed in the previous sections, with 
a straightforward extension to the multivariate case. We evaluate the distribution of K max 



which is the sup of the K' L (a,r) function (see equations [|, [L4|, |15|) using Monte Carlo 
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no detectlor , >' odei. >2crdet. >3crdet. 




0.5 1.0 1.5 2.0 2.5 

Kmax 

FIG. 2: The K max (bivariate) from a Monte Carlo simulation of 500 Gaussian CMB realizations. 
The colored areas indicate the K mSuX intervals where 68%, 95% and 99% of the Gaussian realizations 
are found. If K mSuX is found in one of these intervals, this is taken as a 1, 2 and 3 sigma detection 
of non-Gaussianity respectively. 



simulations of Gaussian CMB realizations. In figures ([!]), (0) and (|3]) we show the results 
for the uni-, bi- and trivariate case. We now define the la, 2a and 3a detection levels as the 
-Kmax values over which we had 68%, 95% and 99% of the hits in the Gaussian simulations 
respectively. These detection levels are shown as shades in the figures. We consider the 
resolution L = 500; to approximate the suprema, we choose a grid of 30 points over the a- 
space, which makes the numerical implementation of our Monte Carlo procedures extremely 
fast. We have some evidence, however, that our results may be to some extent improved for 
finer grids. 

The results are reported in Tables H, [Ij and [TJ. We considered the univariate test 



15 




1.0 1.5 2.0 2.5 

Kmax 

FIG. 3: The -fT max (trivariate) from a Monte Carlo simulation of 100 Gaussian CMB realizations. 
The colored areas indicate the K mSuX intervals where 68%, 95% and 99% of the Gaussian realizations 
are found. If K mSuX is found in one of these intervals, this is taken as a 1, 2 and 3 sigma detection 
of non-Gaussianity respectively. 

described in Section O, and then the multivariate version of Section [TIT]; for the latter, 
we focussed on bivariate and trivariate circumstances, with Al = 250, Am = and 
A£i — Ami — 249, A£ 2 = Am 2 = 250, respectively (see eq.[T^). For this model, the 
performance of the three procedures is certainly very satisfactory. The power of the tests, 
i.e. the percentage of simulations where non-Gaussianity is correctly detected is a monotonic 
function of f3 (the amount of non-Gaussianity) as expected. For instance for the univariate 
test that can consistently detect as little as 6% non-Gaussianity in the map, at the 2a level. 
A map of a model with 6% non-Gaussian distributions are shown in figure (|]) and its his- 
togram in figure (|5]). In pixel space this model seems to be almost identical to a Gaussian 
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model. 

As a second alternative (Model 2), we try to mimic cosmic string models by generat- 
ing 10 4 string-like features with randomly(Gaussian) varying length and temperature and 
superimposing them on a Gaussian background, according to the expression 

T(9, <p) = (1 - P)T G (6, if) + PT S (6, <p), (18) 

where T s (8,(p) is a pure string map, T G (8,(p) is a Gaussian CMB map. The two fields are 
independent. The percentage of non-Gaussianity, /jvc is again defined by (|T7|). We view 
([jjD as a model where non-Gaussianity is strong in pixel space, and we thus expect our 
results not to be as good as for the previous example. An inspection of the Tables |rV|^Vlll| . 
suggests that the performance of the univariate test gives weaker results than what was the 
case for the non-Gaussianity generated in spherical harmonic space: at least 30% of non- 
Gaussianity is needed to ensure a mere 50% detections at la. However, the performance of 
the bivariate and trivariate procedures is more promising, as non-Gaussianity can be detected 
for percentage of non-Gaussianity around 10/15% at the la level. The effect of varying Am 
is noticeable, but not extraordinary; likewise, some unreported simulations for other values 
of A£ (=200,150,100) show similar outcomes (on the other hand, the performance is worse 
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-4 -2 2 4 

u 

FIG. 5: The histogram of the pixels of a Gaussian map with a 5.9% non-Gaussian part, made with 
X 2 distributed ai m . The solid curve is a Gaussian fit to the histogram. 

for very small A£, i.e. 1 to 50). In Figure @ we show a map of a 15.5% strings. The map is 
looking Gaussian except for a few strong out-layers which is seen as a tail in the histogram 
(Figure (0)). This map is at the limit of the 2a detection level (for detection in 50% of the 
cases) for our tests, but it seems that a test in pixel space may here give a stronger detection 
if the out-layers are not misinterpreted as foregrounds or noise. 
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-3.52e-01 ^^^^^^^^^^k 3.52e-01 mK 



FIG. 6: A Gaussian map (monopole and dipole removed) with a 15.5% non-Gaussian part, made of 
string like objects of varying length and temperature. The map has a few pixels which have a much 
higher (positive) temperature than the rest of the map. These strong out-layers are here damped 
so that they don't dominate the map, making us unable to distinguish the other fluctuations. All 
pixels above T max = |T min | are set to T max . 



TABLE I: MODEL 1, UNIVARIATE TEST 



(3 (/ng) 


0.20 (5.9%) 


0.22 (7.3%) 


0.24 (9%) 


0.26 (11%) 


0.28 (13%) 


0.30 (15.5%) 


la 


78% 


96% 


100% 


100% 


100% 


100% 


2a 


57% 


82% 


99% 


100% 


100% 


100% 


3a 


39% 


67% 


95% 


100% 


100% 


100% 



TABLE II: MODEL 1, BIVARIATE TEST 



(3 (/ng) 


0.20 (5.9%) 


0.22 (7.3%) 


0.24 (9%) 


0.26 (11%) 


0.28 (13%) 


0.30 (15.5%) 


la 


61% 


68% 


90% 


99% 


100% 


100% 


2a 


21% 


31% 


47% 


82% 


97% 


100% 


3a 


6% 


9% 


19% 


53% 


90% 


99% 
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FIG. 7: The histogram of the pixels of the map showed in Figure ([]). This is a Gaussian map with 
a 15.5% non-Gaussian part made of string like objects as described in the caption of Figure @. 
The solid curve is a Gaussian fit to the histogram. As small tail is seen due to the high temperature 
strings. 



V. COMMENTS AND CONCLUSIONS 



We believe the approach advocated here may enjoy some advantages over existing meth- 
ods, and we view it as complementary to geometric approaches in pixel space (for instance, 
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TABLE III: MODEL 1, TRIVARIATE TEST 



P (/ng) 


0.20 (5.9%) 


0.22 (7.3%) 


0.24 (9%) 


0.26 (11%) 


0.28 (13%) 


0.30 (15.5%) 


la 


32% 


54% 


76% 


89% 


96% 


100% 


2a 


5% 


17% 


28% 


67% 


84% 


98% 


3a 


1% 


7% 


12% 


33% 


67% 


95% 



TABLE IV: MODEL 2, UNIVARIATE TEST 



P (/ng) 


0.20 (5.9%) 


0.25 (10%) 


0.30 (15.5%) 


0.40 (30.7%) 


0.50 (50%) 


la 


28% 


33% 


27% 


48% 


71% 


2a 


5% 


5% 


3% 


17% 


41% 


3a 


2% 


2% 


0% 


3% 


25% 


TABLE V: MODEL 2, BIVARIATE TEST,A^ = 250, Am = 


P (/ng) 


0.20 (5.9%) 


0.25 (10%) 


0.30 (15.5%) 


0.40 (30.7%) 


0.50 (50%) 


la 


33% 


41% 


55% 


71% 


83% 


2a 


1% 


17% 


22% 


48% 


57% 


3a 


0% 


5% 


9% 


38% 


44% 


TABLE VI: MODEL 2, TRIVARIATE TEST, A^ = 249, At 2 = 250, Ami 


= Am 2 = 


Ung) 


0.20 (5.9%) 


0.25 (10%) 


0.30 (15.5%) 


0.40 (30.7%) 


0.50 (50%) 


la 


30% 


31% 


51% 


72% 


89% 


2a 


6% 


17% 


27% 


56% 


72% 


3a 


1% 


5% 


11% 


46% 


60% 



TABLE VII: MODEL 2, BIVARIATE TEST,A^ = 250, Am = 250 



P (/ng) 


0.20 (5.9%) 


0.25 (10%) 


0.30 (15.5%) 


0.40 (30.7%) 


0.50 (50%) 


1.0 (100%) 


la 


46% 


57% 


73% 


80% 


96% 


100% 


2a 


12% 


27% 


42% 


65% 


84% 


98% 


3a 


4% 


18% 


26% 


57% 


77% 


96% 
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TABLE VIII: MODEL 2, TRIVARIATE TEST, A£ t = 249, A£ 2 = 250, Ami = 249, Am 2 = 250 



P (ha) 


0.20 (5.9%) 


0.25 (10%) 


0.30 (15.5%) 


0.40 (30.7%) 


0.50 (50%) 


1.0 (100%) 


la 


39% 


50% 


65% 


80% 


90% 


100% 


2a 


14% 


27% 


46% 


72% 


86% 


100% 


3a 


6% 


10% 


32% 


61% 


81% 


100% 



methods based an Minkowski functional, local curvature, or other topological properties). 
Our proposal allows for a rigorous asymptotic theory; it is completely model free; it pro- 
vides information not only on the existence of non-Gaussianity, but also on its location in 
the space of multipoles; the effect of estimated parameters is carefully accounted for; given 
that the asymptotic behavior of the field has been thoroughly investigated, many other 
procedures, further than those we considered here (Kolmogorov-Smirnov and Cramer- Von 
Mises tests) can be immediately implemented; the extension to an arbitrary number of 
rows is in principle straightforward (although computationally burdensome), whereas for 
instance the explicit form of higher order cumulant spectra needs to be derived analytically 
in a case by case fashion. Furthermore, our analysis of the distributional properties of the 
spherical harmonic coefficients may have some independent interest in other areas of CMB 
investigation. 

Our approach is clearly related to the analysis of higher-order cumulant spectra (such 
as the bispectrum and trispectrum) in harmonic space. Although the bi- and trispectrum 
have been very widely used in empirical work, their power properties against a variety of 
non-Gaussian models do not seem to have been very much investigated. We conjecture that 
our procedure may enjoy better power properties than higher order spectra in a number of 
circumstances; heuristically, the bispectrum and the trispectrum search for non-Gaussian 
features on the ae m by focusing essentially on their skewness and kurtosis, whereas the 
method we advocate here probe their whole multivariate distribution. The mutual interplay 
between these different methodologies may lead to improvements in both directions: for 
instance, it seems possible to model the bispectrum and trispectrum evaluated at different 
multipoles as processes indexed by some r, much the same way as we did here to combine 
the information from empirical processes G4...4 into a single statistic K^. This may allow 
for a more rigorous analysis in the aggregate, in order to understand whether a single or a 
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few high values are to be considered significant, over a set of L statistics. On the other hand, 
incorporation of some selection rules (e.g., the Wigner's 3j coefficients) into our analysis may 
help us to exploit the isotropic nature of the field to probe non-Gaussianity more efficiently 
In this paper we have tested the method for two non-Gaussian models. One which was 
created in spherical harmonic space and the other which was created in pixel space. For 
the first model, we detect non-Gaussianity at a 2a level (about 50% of the times) with 
the univariate test even when the test model contains only 5.9% non-Gaussianity. In this 
case the map is very similar to a Gaussian map. Our second non-Gaussian test model was 
generated in pixel space. In this case we had 2a detections (about 50% of the times) with 
15.5% non-Gaussianity. The map shows that this kind of non-Gaussianity may be more 
easily detected in pixel space. When taking into account realistic effects like detector noise 
and galactic cuts, the method has so far shown promising results. This will be discussed 



further together with tests on realistic non-Gaussian maps in [31 
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APPENDIX A 



The purpose of this appendix is to derive analytically some of the results concerning 
the bias due to estimated parameters in our procedure; more details can be found in |28f . 
Heuristically, we are going to justify the appearance of 'extra' terms in equations (|), ([13]) 
and (|T5|). To this aim, note first that 

l(u em <a) ee lf%£<-21og(l- a ; 
V Ci 

2\a em \ 2 < 21og(l-a) 



and 



I a>eo 1 2 2 1 an \ z 2 1 ae m \ z \ a „ / 1 





12 


2| 




I 2 


-—£ 




2' '"' W 

Z— /m= 


-—t 




2 / 



D(-, !,...,!), (Al) 



E£ I 12 ' V"^ I 1 2 ' I 12 I ~ Xc 2 

— denoting equality in distribution and 1, 1) a so-called Dirichlet distribution with 



parameters (~,1,...,1) (see Ref.p2|). 

In the sequel, we shall write for brevity 



c \o-m\ 2 c 2|a^ m | 2 

s£o — 737 ; — , Kim 



Et I 12 s^ 1 I 12 

m=— i l a ^"il 2^ik=-e \ a ^m\ 

and 

~ «fe/ ^r 1 ^) «fa/ ®2 l (ai) -21og(l-aj) . 

a * - wr ' a * - wr - — 27Ti — ' 3 ~ 1,2 • 

1. A general result on asymptotic bias 

We select a number p < £ out of £ elements (£a, & P ); of course the univariate, bivariate, 
and trivariate cases follow by setting p equal to 1, 2 and 3, respectively. From PS] we know 
that 

+ D{l,l,...,l-p+~) 
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Note that 

v 

i=l 

= (l{u emi < «i)l(^ m2 < a 2 )...l{ue mp < Oip)) 
= (l(ue mi < a>i,u em2 < a 2 , •••,^m p < a P )) 
= (I (&i < aie, £t P < apt)) 
= P (&i < oiu, ...,£e P < otpt) . 

The technical result we need is then the following (see also f28|j): 
Lemma Al Let T p be the class of all 2 P subsets 7 of {1, 2, ...,p}. For any < au, a p e < 1 
such that Yui=i a iz — ^-> we have 

P(&i<oiu,..,tie P <a pe )= $^(-l) #(7) (l-^a,, 

where #(7) denotes the number of elements of 7. 

Some very simple special cases of Lemma Al, which are of direct interest for the arguments 
below and can be checked by direct evaluation of the corresponding integrals, are as follows: 
forp=l, 7 = {0},{l} 

P(Za<a ie ) = l-(l-a u y- 1 / 2 , 
andforp = 2, 7 = {0},{l},{2},{l,2} 

P (61 < au, in < a 2i ) = 1 - (1 - a u f- 1/2 - (1 - a 2£ f- 1/2 + (1 - a u - a 2e ) £ - 1/2 . 




Proof of Lemma Al We shall give the proof by induction. For p = 1, we have immediately 

/3(a,b) denoting a Beta-distributed random variables with parameters (a, b). Therefore, for 
< a u < 1, 

P < a u ) = r ^^ 1/2) ; : [""(I - u)'- 3 ' 2 du = 1 - (1 - a u y- 1/2 . 

v«i _ u) r(i)r(^-i/2) J K ' 
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For clarity of exposition, we consider explicitly also p = 2, where we obtain, taking into 
accounts that the conditional distribution of — is (3(1, £ — 3/2), 



]_ ran 

P (&i < au, £fi < a n ) = (£ - -) jf (1 - a:] 



£-3/2 



1-(1 



ft 2£ ^-3/2 



1 — X' 



dx 



an 



{£--) / (l-x) e - 3/2 dx -(£--) / (1 -x-a 2e ) e ~ 3/2 d 



"if 



/./• 



= 1 - (1 - au) 1 - 11 * " (1 " M £ ~ 1/2 + (1 " «ii - « 2 ,r i/2 
For general p > 2, we have (see Ref. f32H) 

P{ili < a u , ...,£t p < 
T(£+l/2) [ au 



r(i- P +i/2) y 



e-p-1/2 



i=l 



du\...du p , 



which becomes, after the change of variables Wj = 1^/(1 — w p ), i = 1, 2, — 1 

p-1 \ «-p+V2 



r(£ + i/2) p** 



r(£- P + i/2) J 

T(£+l/2) [°* 



(i - Up y--^ 



a ie /(l-u 3 ) pa p _ u /(l-u 3 ) 



t=l 



dw\...dw v 



T{£-l/2) 



l ap \i- Up y- 3/2 p(ta< T 



au 



< 



u 3 



a P i 
1 - u 3 



du?, 



r(£ + i/2) 

T{£- 1/2) 



3j 



(A2) 



the last equality following from the induction step. Now ( |A2| ) becomes 



£-3/2 



M 3 - E a j£ 



£-1/2 



du 3 

£-1/2 ' 



7'6r p _i 



J67 



1 - a pi - 22aj£ 

£7' 



£(-!)*<■» 1-J> 



£-1/2 



JG7 



as 



required, by posing 7 = (7' U {p}), and because T p = T p _i U {(7',^) : 7' e r p _i} . 



□ 



From Lemma Al, and for £ suitably large 



2£+l J' ^ 2£+l 
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and it can be readily checked that 



Finally, note that, as L — > oo 

[Lr] 



1 f 1 1 

> ; ~ b(a)\/r / —=du = 2b{a)y/r . 



The calculations for the multidimensional cases follow in very much the same way. 



2. The bias for the m = term 

We show now that the bias for the term corresponding to m = is asymptotically 
negligible. Note first that — and define 

El=l l a *n| 2 

which has a Snedecor F-distribution, with 1 and 21 degrees of freedom. We have 

\a m \ 2 Tt 



whence, after straightforward calculations 

; , |a £ 2 ■ !>-" (A3) 
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l«£0| 2 + 2Em=l l^m| 

P | * i ( ^ ) <a)-a (A4) 

V ~ 2£+l-^\a) J V ; 

< 2£ ^ 1 ff (a) ) - p ( T * < $r») + < $ rV)) - <* 

P < TTTT^St^v^ - P ^ < (A6) 



l + ( l_$ r i (a))/2 £ 
+P(T e <®i\a))-a. (A7) 

The first component fl^|) for £ large enough is bounded by 

PI E>) <T,< ^ 



1 + |1 - ^\a)\/2l ~ ~ 1 - |1 - $r 1 («)|/2^ 
P v 1 \ J < JT e < v 1 V/ A8 

0(7). 
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because \fTt has a Student's t-distribution, the latter has a bounded density function and 
the length of the interval in ( |A8|) is of order £~ l . For ( |A7|) , we have that 



P(T i <<$>-\a))-a = 2^ 



f t -i(u)du 



(A9) 



fi(x) and ip(x) denoting, respectively, the density functions of a Student's t distribution 
with 2£ degrees of freedom, i.e. 

r((f + i)/2) 1 1 
Jey > T{£/2) ^fa(i + u 2 /ty e + i y 2 ' 

and a standard Gaussian density. Thus ( |A9[ ) is bounded by 

r((£ + i)/2) i i 



+2 



e/2T(£/2) V2n 

r((£ + i)/2) i 



2n 



exp(— u /2)du 



{(1 + u 2 /£)- £/2 - exp(-u 2 /2)}du . 



For ( gTL 



we use 



sup 

0<u<x 



1 + 



■it 



2\ -(«+l)/2 



exp(-w 2 /2) 



0( 



and finally, (|A10|) follows from the well-known property 

r ^_ +1 ^-lUo(l),as/^oo. 
/£/2T(£/2) J r ' 

Thus it follows that (|A10|) is 0(l/£) as well, and the proof is completed. 



(A10) 
(All) 
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